function Rho = denfunc
N = 128;
ntotal = N^2;
epsilon = 1/N;
lc = 8;

h = (1/N)/epsilon;
[x,y] = ndgrid( (0:(N-1))/N, (0:(N-1))/N );
% V = 1d-2*sin(2*pi / lc * (x/epsilon)) .* ...
  % sin(2*pi / lc * (y/epsilon) ) + 1d-3 * (rand(size(x))-0.5) .* ...
  % (rand(size(y))-0.5) ;
% V = zeros( N, N );
V = 1d-1 * reshape(rand(ntotal,1), N, N);
nx = N;
ny = N;
display('Generate a 2-d tight binding matrix');
H = sparse(nx*nx, nx*nx);
for i1 = 1 : nx
  for j1 = 1 : ny
    pos = nx*(i1-1) + j1;
    pos1 = nx * (i1-1) + mod(j1,ny)+1;
    pos2 = nx * (i1-1) + mod(j1-2,ny)+1;
    pos3 = nx * (mod(i1,nx)) + j1;
    pos4 = nx * (mod(i1-2,nx)) + j1;
    H(pos, pos1) = -1;
    H(pos, pos2) = -1;
    H(pos, pos3) = -1;
    H(pos, pos4) = -1;
    H(pos, pos) = 4;
  end
end
H = H / h^2;
for j1 = 1 : nx
 for i1 = 1 : ny
   pos = nx*(j1-1) + i1;
   H(pos, pos) = H(pos, pos) + V(j1,i1);
 end
end

display('Diagonalizing');
tic
NSite = (N/lc)^2;
opts.disp = 0;
[EF, D] = eigs( H, NSite+10, 'sm', opts );
toc
D = diag(D);
SD = sort(D);


beta = 500;
mu = SD(64);
% mu = SD(16);
Rho = real(2d0 * abs(EF).^2 * ( 1./ (exp(beta*(D-mu))+1d0) ) );
Rho = reshape( Rho, N, N );
